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ABSTRACT 



Aims. The time evolution of the electromagnetic emission from blazars, in particular high frequency peaked sources (HBLs), displays 
irregular activity not yet understood. In this work we report a methodology capable of characterizing the time behavior of these vari- 
able objects. 

Methods. The Maximum Likelihood Blocks (MLBs) is a model-independent estimator which sub-divides the light curve into time 
blocks, whose length and amplitude are compatible with states of constant emission rate of the observed source. The MLBs yields the 
statistical significance in the rate variations and strongly suppresses the noise fluctuations in the light curves. 

Results. We apply the MLBs for the first time on the long term X-ray light curves (RXTE/ASM) of Mkn 421, Mkn 501, lES 1959+650 
and lES 2155-304, which consist of more than 10 years of observational data (1996-2007). Using the MLBs interpretation of 
RXTE/ASM data, the integrated time flux distribution is determined for each single source considered. We identify in these distribu- 
tions the characteristic level as well as the flaring states of the blazars All the distributions show a significant component at negative 
flux values, most probably caused by an uncertainty in the background subtraction and by intrinsic fluctuations of RXTE/ASM. This 
effect interests in particular short time observations. In order to quantify the probability that the intrinsic fluctuations give rise to a 
false identification of a flare, we study a population of very faint sources and their integrated time flux distribution. We determine 
duty cycle or fraction of time a source spent in the flaring state of the source Mkn 421, Mkn 501, lES 1959+650 and lES 2155-304. 
Moreover, we study the random coincidences between flares and generic sporadic events such as high energy neutrinos or flares in 
other wavelengths. 
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1. Introduction 

Blazars are defined as Active Galactic Nuclei (AGNs) domi- 
nated by a highly variable component of non-thermal radiation 
produced in relativistic jets pointed close to the line of sight 
(iBegelman et al. 198 41: lUrr y et al.l995 ). One of their main char- 
acteristics is the flux variability on different time-scales: from 
fast flares lasting few minutes to high states of several months. 
Blazars are considered to be sites of energetic particle produc- 
tion and potential sources of cosmic rays up to energies of at 
least 10'^ eV. 

The standard blazar spectral energy distribution (SED) shows 
two prevalent components: a hump at low-energy which peaks 
in the frequency range between infrared and X-ray bands, and 
a second hump at higher energy, proportionally shifted in the 
range from MeV up to TeV y-rays. Two potential scenar- 
ios, the so-called leptonic and hadronic ones, have been pro- 
posed in order to model the SED. In leptonic mod els (e. g 
iJones et al.l974: GhiselHni et al.l996; Mastichiadis et al. 19971) . 
synchrotron emission from relativistic electrons is responsi- 
ble for the first hump. Electrons in the jet plasma up-scatter 
low-energy photons to high energies via Inverse Compton, 
producing the second hump. In this scheme, the same 
electron pop ulation produces both components. In hadronic 
models (e.g Mannheim and Biermann 19891; iMannheim 19931: 
lAharonian 2000h protons are accelerated in the jet together 



with electrons. The synchrotron radiation produced by pri- 
mary and proton-induced electrons contribute to the low-energy 
component. High-energy radiation originates from photo-meson 
interactions and from proton and muon synchrotron radia- 
tion. A comprehensive description of a Monte Carlo simu- 
lation of a stationary synchrotron proton blazar model, in- 
cluding all relevant em ission processes, can be found in 
jMiicke and Prothereo 2000,) . In hadronic models, y-ray produc- 
tion by pion photo-production results in simultaneous neutrino 
production. The decay of charged pions is the main neutrino pro- 
duction channel as discussed in (Miicke et al. 20031). 
The detection of very high energy neutrinos coming from blazars 
would be an unambiguous proof of the existence of baryonic 
loaded outflows and would indicate that blaz ars accelerate high 
energy cosmi c rays. Neutrino telescop es (e.g. A hrens et al.2004t 
Ackermann et al. 20071; lAntipin et all007) until now have not 
detected any extraterrestrial sour ce of neutrinos in t he TeV-PeV 
energy region. As discussed in jRachen et al.I998l) . "the tran- 
sience of energetic emission could improve the association of de- 
tected neutrinos with their putative sources, because one could 
use both arrival direction and arrival time information, allowing 
statistically significant statements even for total fluxes below the 
background level". This is true under the assumption that neu- 
trino production in HBLs is subjected to the same mechanisms at 
the base of the electromagnetic activity. Consequently, neutrino 
production and electromagnetic activity should show the same 
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time modulation. The observation of time coincidences between 
electromagnetic flares and rare events, like neutrinos, represents 
a natural test to the hadronic scenario. 

The main requirement to this approach is a clear definition and 
classification of the states of activity of the observed source. 
In this paper we discuss a procedure able to identify char- 
acteristic and flare states in a light curve. The estimator that 
best fits our requirements is the Maximum Likelihood Blocks 
(MLBs), since it is model-independent, it has been designed 
to identify blocks of data with a constant rate in variable pe- 
riods, and it provides a statistical significance for each block. 
To test our approach, we perform a complete and detailed 
analysis on RXTE/ASM X-ray light curves. In particular, we 
analyze data from the brightest High energy peaked BLLacs 
(HBLs) (iGiommi and Padovani 19941) : Mkn 421, Mkn 501, 
lES 1959-H650 and lES 2155-304. In the first part of this pa- 
per we describe the MLBs and how to separate flares from the 
characteristic level. Moreover, we introduce a definition for the 
duty cycle of the source. In the second part, we discuss the ap- 
plication of the method on RXTE/ASM data. 

2. Methods 

A variety of methods are used in astrophysics in order to as- 
sess the variability of a source and to qualify the character of 
the variability (periodic, correlated etc). It is not our intention 
here to review these methods. Each method is designed for a 
specific purpose. Often, data are affected by large uncertainties 
or the data spacing is rather inhomogeneous. The driving fac- 
tors for the selection of a method are the goals of the analysis 
and the quality of the data to be analyzed. In our case we need a 
method that addresses the variability issue on light curves which 
are unevenly spaced and have short and long breaks, takes into 
consideration the statistical errors and possible unknown instru- 
mental fluctuations on the measurements and gives a represen- 
tation of the light curve in term of periods in which the data 
points are compatible with a constant level. A method that could 
satisfy these requirements is the Maximum Likelihood Blocks. 
The entire data analysis rep orted here is performed in ROOT 
dBrun and Rademakers 1997 ), an object-oriented data analysis 
framework. The only exception is for the Maximum Likelihood 
Blocks algorithm which is currently an IDL based program. 

2.1. Representation of the light curve: Maximum Likelihood 
Blocks 

The methods used in the study of temporal variability depend 
strongly on the nature of the available data and of the signal of 
interest. In all cases, the most basic step is the classification of 
the time-series as "constant" or "variable". Suitable and widely 
used statistical tests include the Kolmogorov-Smirnov test for 
time-tagged data (e.g. the arrival times of X- and y-ray pho- 
tons) and the x'^^ test for binned data (e.g. binned photon ar- 
rival times or optical magnitudes). The next step in the analy- 
sis of light curves is the characterization of their "shape". We 
will use a simple and model-independent approach that aims at 
dividing the light curve into time intervals in which the source 
emission is compatible with a constant level. An algorithm based 
on Bayesian statistics that performs such a se gmentation fo r 
data of diff'erent natures was presented by (S98, IScargle 1998]) . 
In its form for time-tagged data, this algorithm was used for 
example to characterize the X-ray light curves of young pre- 
main sequence stars ob served by the Chandra Orion Ultra-Deep 
Project (COUP, .Getman et al.2005h in the Orion Nebula Cluster 



(ONC). A modified version of the S98 algorithm, based on a 
Maximum Likelihood rather than a Bayesian approach, was re- 
cently employed in other studies of stellar X-ray lig ht curves 
dWolk et al.20051: iFavata et al.2005l: [Stelzer et al.2006l) . We will 
refer to this algorithm as the Maximum Likelihood Blocks 
(MLBs) and we introduce here a variant that is suitable for the 
analysis of binned light curves. 

Our algorithm is derived from the one presented by S98. Like 
S98, we tackle the problem of finding the best piecewise rep- 
resentation of a binned light curve in an iterative (and approxi- 
mate) way: we begin by testing the data against a constant-flux 
model. If the model does not represent the data adequately we 
split the light curve into two parts at the most likely "change 
point". We then repeat these two steps recursively on the re- 
sulting segments until all segments are compatible with constant 
emission. The fundamental difference with the algorithms pre- 
sented by S98 lies in the statistics used to test if a light curve 
is variable and to find the most likely change point: rather than 
"marginal likelihoods" and "Bayes factors" (e.g. Eq. 7 and 48 in 
S98) we employ simple Likelihood functions, i.e. the probability 
densities of obtaining the observed data set given a parametric 
model. 

As mentioned above the algorithm was first applied to time- 
tagged data. It is here adapted to the different statistical prop- 
erties of binned time series. Our lightcurves can be described 
as a series of independent flux measurements, r,, each normally 
distributed about their mean values with standard deviations cr, . 
The likelihood of a parametrized model, M, of the lightcurve is 
maximized by minimizing the ^2 - J](r; - r/[M])2/cr/2), where 
r,[M] are the model-predicted fluxes. In our case the model M 
is either the single-segment representation or one of the possible 
two-segment representations of the light curve. We will refer to 
these models, respectively, as "1", specified by one parameter, 
the constant flux level, and "2(j)" specified by three parameters: 
the change point "(j)" (more specifically the index of the last 
point in the first of the two segments) and the two flux levels. In 
this notation our algorithm reduces to: /) splitting the light curve 
if the minimum ;t'i2 is such that the probability of obtaining a 
larger value is lower than a user defined confidence threshold 
(e.g. 1% or 0.1%); ii) choosing the change point, jcp, as the one 
that minimizes 

2.2. Interpretation of the light curve: flares versus 
characteristic level 

The goals of this analysis are the identification of the vari- 
ous levels of activity of a source and the separation between 
bursting events (flares) and steady state period(s) (characteris- 
tic level(s)). Sometimes periods of no variable activity are de- 
fin ed in the literatur e as "quiescent". As discussed for example 
m ( Wolk et al.20"05 l). the meaning of quiescent emission is am- 
biguous. An apparently quiescent level can be due to a super- 
position of numerous unresolved flares. Quiescent, as defined as 
inactive, is therefore not appropriate to describe the level of ac- 
tivity in which the source spends most of the time. We define the 
characteristic level as Rci,ar and the spread around it cTchar- In 
order to determine the value of Rdwr we construct the distribu- 
tion of the amplitude r, and the duration of the single block T,. 
We call this integrated time(T)-flux( r) distribution based on the 
MLBs interpretation (B): Tsir). This provides the distribution of 
the total amount of time the source passes in a particular activity 



' The minimization with respect to the flux levels is trivial and re- 
duces to choosing the mean flux in the given time interval. 
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state. The threshold above which a flux state is defined sls flare 
is then defined as R^cr = {Rchar + NcTchar)- Depending on A^, the 
probability that a selected flare state is caused by a fluctuation 
of the characteristic level, by an instrumental fluctuation or by a 
real enhancement of the photon emission from the source can be 
fully assessed. 

On the base of this definition of flares we can determine as well 
the frequency of flare states or duty cycle Df^u- as: 



Durr — 



r TB(r)dr 
rmr)dr 



(1) 



In Sect. 4, the application of this method to RXTE/ASM data 
are reported. 

3. Data 

The All-Sky Monitor (ASM) on board of the Rossi X-ray Timing 
Explorer (RXTE) has been monitoring the X-ray sky routinely 
since March 1996. During each orbit up to 80% of the sky is sur- 
veyed to a depth of 20-100 mcrab. A source is observed roughly 
10 times a day. A set of linear least-square fits over 90 sec- 
onds observation periods, by one of the three Scanning Shadow 
Cameras, yields the source intensities in four energy bands (1.5- 

3, 3-5, 5-12, and 1.5-12 keV). The intensities are usually given 
in units of the count rates expected if the sources were at the 
center of the field of view in one of the cameras. In 1.5-12 
keV band, the Crab Nebula flux corresponds to about 75 ASM 
co unts per second. A detailed description of ASM can be found 
in jLevine et al. 19961) . RXTE standard data products are col- 
lected directly from the HEASARC database. 

We concentrate our study on RXTE/ASM data because this pro- 
vides the longest light curves in X-ray of Mkn 421, Mkn 501, 
lES 1959-H650 and lES 2155-304. However, for these kind of 
sources, the RXTE/ASM sensitivity is limited and data are af- 
fected by large errors. Moreover, the resolution and the back- 
ground level of ASM observations depend on the Sun contami- 
nation or back-scattered solar X-rays and on the dete ctor stabil- 
ity along the 10 years of data taking, dWen et al.2006l) . 

4. Results 

The results of the application of the MLBs to the RXTE/ASM 
data for the four HBLs considered are reported in Fig. [1] and 
in Fig. |2] Each change point identified by the algorithm has a 
statistical significance of at least 3cr. This operation is still not 
enough in order to characterize the behavior of a source. We 
haven't yet quantified the characteristic noise of RXTE/ASM, 
then we can not distinguish if the change points identified by 
the MLBs are of instrumental nature or of a physics nature. In 
order to distinguish between these two scenarios and study the 
behavior of sources, we construct the integrated time(T)-flux(r) 
distribution Tgir), as discussed in Sect. 2.2. We first construct 
the Tgir) for a set of very faint sources in order to study the 
intrinsic fluctuations of the instrument and finally we determine 
the rB(r) for Mkn 421, Mkn 501, lES 1959+650 and lES 2155- 
304. On the base of these distributions we identify the states that 
can be considered flares with good confidence. 

4. 1. RXTE/ASM intrinsic fluctuations 

In Fig. [T] we notice that the MLBs identify not only signifi- 
cant change points at positive amplitudes but also at negative 



ones. These negative fluctuations can be caused by uncertainties 
in the background subtraction and by intrinsic fluctuations of 
RXTE/ASM. In order to characterize such a component and its 
effect on the definition of flares, we have analyzed RXTE/ASM 
light curves for a set of very faint sources, since these are ex- 
pected to spend most of their time at a flux level well below 
the ASM sensitivity. The sources considered, reported in Tab.[Tl 
have a low X-ray monochromatic average emission (less then 
0.6 yuJy at IkeV), are randomly distributed in the sky and are 
at various redshifts. The flux distributions of the faint sources 
are all normal distributions, as expected for a random instru- 
mental noise. On average, the normal distributions peak at rate 
r w 0.1 ASM c/s and have a standard deviation of r ^ 0.3 ASM 
c/s. Since the studied faint sources show similar flux distribu- 
tions, we will use in the next just one of them for compari- 
son; the source used is PKS 0118-272 and represent the aver- 
age faint source in our sample. In Fig. |3] the flux distribution 
of PKS 0118-272 is compared with the flux distribution of the 
HBLs considered in this work. The distributions are normalized 
using the areas under the negative flux tails. In this way, we es- 
timate the fraction of the HBLs flux distributions caused by the 
intrinsic fluctuations of RXTE/ASM (S no-)- 

4.2. RXTE/ASM flare states 

The Teir) for the Mkn 421, Mkn 501, lES 1959+650 and 
lES 2155-304 are reported in Fig. [3] All the distributions 
differ significantly from a normal distribution indicating that 
RXTE/ASM is indeed sensitive to high activity states of the 
sources considered. For all the HBLs considered we observe 
that the flux distribution shows a peak above the pure back- 
ground distribution. We define the central peak value (that cor- 
responds also to the mode of the flux distribution) and its stan- 
dard deviation as Rdwr and cTchar- More sophisticated fitting pro- 
cedures have been applied but, given the quality of the data 
of RXTE/ASM, they did not provide a more precise estima- 
tion of Rchar and cTchar- We observe that Rchar is next to the 
detection threshold of RXTE/ASM. As discussed in Sect. 2.2, 
Rnct = (Rchar + NcTchar) IS the threshold above which a flux state 
is considered a flare. 

Using this definition of flares we determine the frequency of flare 
states or duty cycle D^cr as described in Sect. 2.2. In Tab.|2]we 
report D^cr for the HBLs considered and the cases of N-l and 
N-3. For the specific case of RXTE/ASM, we calculate as well 
the intrinsic fluctuation S no- that affects the duty cycle as: 



Sno- - 



(2) 



where Tb^^J/) is the flux distribution of a faint source (in this 
work PKS 01 18-272). Results are reported in Tab.|2] 



4.3. Example of application: correlation study between 
electromagnetic flares and neutrinos 

As anticipated in the introduction, the study of the physics 
of HBLs develops through different approaches. One of the 
most frequently used sees the study of flare correlation among 
different waveleng ths, for example X-ray and TeV-y rays 
(iMaraschi et al.l99 9). A study of the correlation among differ- 
ent mess engers such as photon s and neutrinos is a more recent 
interest ( Achterberg et al.2005l) and has been the motivation of 
this work. The significance of such correlations can be assessed 
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only when the frequency of the electromagnetic flare states is 
determined, for example following the procedure described in 
this paper. In order to illustrate such a case, we study the dis- 
tribution of coincidences between RXTE/ASM flare states and 
a set of N neutrinos or equivalent sporadic events. The flares 
are selected following the procedure described above for 3crchar 
and the N neutrino events are uniformly distributed in the entire 
time period considered in this paper, approximately 10 years. 
The distributions of coincidences between the RXTE/ASM flare 
periods for Mkn 421, Mkn 501 and IBS 1959+650 and the neu- 
trino events are reported in Fig.|4] Depending on the number N 
of sporadic events we are able to determine the number of ran- 
dom coincidences multiplying the probability of a flare, or duty 
cycle of the source, and N. Once the random coincidence among 
flares and neutrinos is determined, then the statistical meaning of 
experimentally observed coincidences between flares and neutri- 
nos can be determined. This eventually can hint at the associa- 
tion of detected neutrinos with an astronomical source even if N 
is below or at the level of the expected background. 

5. Conclusions 

The X-ray time behavior of Mkn 421, lES 1959+650, Mkn 501 
and lES 2155-304 has been characterized using approximately 
10 years of data from RXTE/ASM. The characteristic level and 
flaring states have been defined and values are reported in Tab.|2l 
Mkn 421 is the source that flares most often amongst those stud- 
ied: for ~ 40% of the RXTE/ASM observations the source was 
in an active state (above Rur) and for ~ 18% of these it was in 
a very high state (above R^a-)- The probability that a flare state 
is caused by a fluctuation of the instrumental noise is marginal. 
This confirms the well known fact that Mkn 421 is an extremely 
variable HBL and quantifies for the first time its duty cycle even 
if this is only valid for RXTE/ASM. Mkn 501 flares often at low 
rates (~ 26%) and less often at high rates (~ 10%). Also in this 
case, the intrinsic fluctuations do not significantly affect the flare 
states. lES 1959+650 flares more rarely in particular at high 
rates (only 2.6%). The systematic component affects one tenth 
of flares. In the case of lES 2155-304 nearly one fourth of flare 
is caused by an intrinsic fluctuation of RXTE/ASM. This simply 
means that this source is at the threshold limit for RXTE/ASM. 
The study of the significance of correlations of flares among dif- 
ferent wavebands and different messengers is foreseen for a fu- 
ture work. 
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Table 1. HBLs and X-ray faint sources considered in this work. The first column contains the X-ray monochromatic emission, the 
second and the third contain the mean and standard deviation of the Gaussian that fit at the best the flux distribution. Last column 
contains the red-shift. The mean value of the RXTE/ASM count rate is reported for the faint sources only since the HBLs are too 
variable for a definition of a mean flux value. For comparison, the resulting Rdwr and <Jchar are reported in Tab.|2] 



Source Name 


X-ray at IkeV (yjjy) 


Mean (ASM c/s) 


Sigma (ASM c/s) 


z 


Mkn421 


10.1 






0.031 


Mkn501 


9.5 






0.034 


lES 1959+650 


9.4 






0.048 


lES 2155-304 


15.2 






0.117 


PKS 0118-272 


0.22 


0.03 


0.25 


0.557 


lES 0235+164 


0.28 


0.09 


0.24 


0.94 


S5 0454+844 


0.03 


0.07 


0.21 


0.112 


PKS 0735+178 


0.20 


0.11 


0.31 


0.424 


PKS 0829+046 


0.08 


0.07 


0.25 


0.174 


S4 0954+65 


0.56 


0.09 


0.21 


0.368 


B2 1147+245 


0.04 


0.10 


0.24 


-0.2 



Table 2. For each HBLs considered we give the characteristic level Rdwr and its standard deviation cTchar- The duty cycle for active 
state Dio- and for very high states D30- is also reported together with the relative intrinsic fluctuations still present above the Rchar+^cr 
or Rchar + 3cr threshold due to RXTE/ASM {S xcr, S 3a-)- 



Source Name 


R,,,ar (ASM c/s) 


CTrfM,- (ASM c/s) 


Dur (%) 


5i.(%) 


£>3. (%) 


(%) 


Mkn421 


0.5 


0.4 


40.3 


1.0 


18.1 


0.5 


Mkn501 


0.4 


0.2 


25.8 


3.0 


9.7 


2.3 


lES 1959+650 


0.3 


0.2 


14.2 


7.1 


2.6 


9.5 


lES 2155-304 


0.2 


0.4 


27.0 


22.1 


4.9 


18.8 
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Time (MJD) 

Fig. 1. Comparison between RXTE/ASM light curve (black points) and the MLBs (red Hne) for Mkn 421, Mkn 501, lES 1959+650, 
lES 2155-304. The period of time reported in the figure is a sub-period respect the one analyzed; data considered in the paper have 
been collected for the period MJD 50200-53698 that corresponds to nearly 10 years. 
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Fig.3. Integrated time flux distribution for Mkn 421 (top-left), Mkn 501 (top-right), lES 1959+650 (bottom-left) and lES 2155-304 
(bottom-right). The distributions are compared with the integrated time flux distribution of a very faint source (dark background). 



